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Abstract We present a visco-elastic coupling model 
between caked spheres, suitable for DEM simulations, 
which incorporates the different loading mechanisms 
(tension, shear, bending, torsion) in a combined manner 
and allows for a derivation of elastic and failure prop- 
erties on a common basis. In pull, shear, and torsion 
failure tests with agglomerates of up to 10000 particles, 
we compare the failure criterion to different approxima- 
tive variants of it, with respect to accuracy and com- 
putational cost. The failure of the agglomerates, which 
behave according to elastic parameters derived from the 
contact elasticity, gives also insight into the relative rel- 
evance of the different load modes. 
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1 Introduction 



X 



cemented contacts 



The research of granular media can be separated into 
two rather distinct fields: There are systems with purely 
repulsive granular interaction ( "dry granular media" , 
e.g. [5]) and media with non- negligible additional at- 
tractions ("wet granular media", e.g. (Oj). While for the 
former the further distinction between granular gases 
(e.g. [TT]) and dense systems can be made, the latter 
is more adequately separated into interactions of per- 
manent or recurring nature (like e.g. van der Waals at- 
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traction or liquid menisci) and those which vanish irre- 
versibly when being overcome (like sinter necks or other 
solid bridges) . In this work, we focus solely on this last 
type. 

Such cemented but breakable contacts form (delib- 
erately or accidentally) in agglomerates of particles. 
The failure behavior of these agglomerates is of prac- 
tical interest in both cases, but also from the point of 
view of Distinct Element Methods (DEMs), modelling 
the hindering of all six relative degrees of freedom of 
two particles in contact is much less studied than the 
well established frictional contact models like the linear 
spring-dash pot, the Hertz contact and its extensions 

(cf. e.g. ^m)- 

For modelling purposes two aspects of the contacts 
have to be taken into account: its elastic response to 
deformations and the maximal load it can bear. The 
elastic properties are usually expressed as set of (often 
decoupled) springs associated to each relative degree 
of motion (cf. e.g. [4l[nil5)). Also for the maximal load 
it is a common practice to employ separate thresholds. 
Exceptions (in the context of granular assemblies) to an 
independent treatment of just tensile and shear load are 
[3jJ:, 15^ who employ empirical or plausible combination 
functions. 

It is thus legitimate to use a simple but self consis- 
tent method to model the cemented contact by a flat, 
elastic cylinder subjected to the Tresca failure criterion 
(section [5]) . The full Tresca criterion is only solvable 
with iterations therefore we compare the full criterion 
to its two simplified variants in DEM simulations of 
the quasi-brittle failure of a cylinder made up of 200 to 
10000 particles (section|31). We conclude with discussing 
the relative relevance of the different loading modes in 
section 14.41 
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2 Modelling breakable cemented contacts 



2.1 Contact model 



The principle modelling idea for the cemented contact 
between two particles is the dominance of the bridge 
with respect to the elastic and failure properties. Hence, 
we envisage the bridge as a flat cylinder (height d smaller 
than diameter 2a) while treating the rest of the two par- 
ticles as solid bodies. 

We set up a local coordinate frame, where the bot- 
tom "plate" of the cylinder is fixed in the x/y- plane 
with its center in the origin, while the upper one is lo- 
cated (in the unstrained state) parallel at z = d (i.e. 
the original contact normal is in the z direction). We 
call the torsion about the z-axis "twist" with angle 4>, 
while we choose the a;-axis of the local frame as the axis 
of the other rotational load, the "tilt" 9 of the upper 
plate. These rotations (being, due to the brittle limit, 
small as to allow for commutation) are accompanied by 
the translational displacements, the pull/push Z and 
the shear X, Y. In total, we describe the load geomet- 
rically by these five parameters. This is independent of 
the actual representation of particle orientation e.g. in 
a DEM code. 

Next, we need to know the elastic response in terms 
of pull/push force Fz, shear forces Fx,Fy, and torques 
TzjTx opposing the twist and tilt, respectively. We avoid 
introducing five arbitrary spring constants, and aim to 
derive them from the elastic behavior of the cylinder, 
instead. 

The stress inside a deformed cylinder (possessing 
elastic modulus E and Poisson ratio v) is in principle 
a text-book problem (cf. e.g. [13]). We directly provide 
the displacement field associated to each loading mode: 



Upush/puii(a;,y, z) 

Ubond(a;,2/, z) 
UshGar(a;,y,z) 



Eira^ 
4r,(l + z/) 



{—VXex — VXGy + ZGz 

{-yze^ + xzey) 



Eira^ 

-^-^ {-2iyxyea; + 2vyze 

-{viy"^ -x^) + z^)ey) 
2 + 2v 



(1) 

(2) 

(3) 
(4) 



From these fields, we can immediately read off the de- 
formation geometry as 



X = e, • Ushcar(0, 0, d) = '^^^^^Fx 

Y = ey U,hcar(0, 0, d) = ^^^-^^Fj, 
Z = Gz ■ Upush/pull(0, 0, d) = — 



e.,, • uto 



Ena'^' 

I Ml + z^ 

(a, 0, d)/a = — — ; — 



-dz{ey ■ Ubond(0,0,d)) 



Ena"^ 

Ad 
E^ 



Tx 



(5) 
(6) 
(7) 
(8) 
(9) 



Summing up (P) to ([U, calculating the correspond- 
ing strain and applying Hook's law yields a stress tensor 
of 



^xy ^xy ^yy " ^ 

_ a^Fx ~ 2Tzy 
"""" ~ 7ra4 

a'^F, + 2T,x 



a^Fx + AT.y 



(10) 
(11) 
(12) 

(13) 



which obeys vanishing divergence and vanishing normal 
stress at the free surface. Integration over the lateral 



surfaces confirms the force F 



-F^e^-f i^^e^ and 



reveals a torque (width respect to the point ezd/2) of 
T = {Tx+Fxd/2)ex+Fyd/2ey+Tzez on the bottom one 
(and thus — F and — T on the top one, of course). The 
additional torque contributions are indeed necessary, 
since we defined the shear mode as keeping bottom and 
top plate parallel. 

With the definition of the spring constant fc„ = 
Et:o? /d (which in principle can be measured in exper- 
iments with two particles |6j), the relation between de- 
formation and reaction force finally reads 



/ Fx \ 




Fy 




Fz 
Tx/a 




Ty/a 




\T./aJ 
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1 

oi 







\ 









4+4i7/ 



Y 

z 

a6 
WJ 



(14) 



When considering the cylinder (the geometry of which 
does not explicitly enter into the contact geometry any- 
way) as being very flat (i.e. d/a — >■ 0), the stiffness ma- 
trix looses the Tj^-line and becomes diagonal. This de- 
coupling is very convenient for DEM modelling. More- 
over, every force/torque component can be supplemented 
by a viscous damping involving the time derivative of 
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the deformation and with coefficients for critical damp- 
ing: 



-F^ = + ktX 



'Fy + ktY 



-Fz = 2yJknmZ + kn Z 
-Tx = a\/ knix 6 + ^'"'^ 



4 



(15) 
(16) 
(17) 

(18) 
(19) 



Here, kt = kn/{2 + 2v) while m, I^^ and Iz are the re- 
duced mass and moments of inertia of the two-particle 
system, respectively. Critical damping is used for con- 
venience when the deformation rate is anyway slow 
enough as to neglect velocity effects. 



2.2 Failure criterion 

Having worked out the elastic response of the bridge, 
we now turn to its failure. The choice of failure criteria 
for condensed matter to be found in the literature is 
vast. We will use the well known Tresca criterionlM, of 
maximal shear stress. From the stress tensor pil)) - (|13l) 
the principle stresses are readily found to be 



7ra^CT2 = 



4 

7ra CT3 



= a^Fzl2 + 2TxV + 



(20) 
(21) 

(22) 



where 



2a\FzTxy + 2Tz{FyX - Fxy)) 



(23) 



Since ai < a2 < crJ3, the maximal shear stress (with 
respect to orientation) is 



0'Trosca(2;, V) 



CT3 - (Ti 



(24) 



and has yet to be maximized with respect to < 
a} and then compared to a material dependent critical 
stress (7*. That means, we are left with a maximization 
of Introducing forces = T^/a, = T^/a and 
dimensionless position x = x/a, y = y/a yields 



§-^=f?ix+^ 
2T, 



4(T2-fT|) 



(25) 



^ We use negative sign for compressive stress. 



where <P^ and 43^ differ only by an irrelevant x/y- 
independent term. Obviously, 3^ is just a quadratic 
form in x and y the center and steepness of which de- 
pend on the load. Unfortunately, when maximizing 3^ 
subjected to the constraint x'^ + y'^ < 1 (which, due 
to the absence of local maxima, can be replaced by 
x'^ + y'^ = 1) one is faced with a quartic, the radicals of 
which exceed any sensible formula length. 

For several special cases, the calculation can be per- 
formed, though: In the simplest case of zero torques, 
CTrcsca bccomes independent of x and y and thus 



na (TTrcsca,! 



= F 



Tresca 



= ^ FS 



F2 

y 



F^/4 (26) 



and can be compared directly to the threshold F^ = 
Tracer*- A less simple case is Fj^ = which lets 3^ take 
on its maximum a,t y — sgn{FzTx — 2FxTz), producing 
the value 



, , F2 \FzTx-2FxTz\ T^ + T^ 



(27) 



Yet another illuminating combination, namely Fx 
Fz = 0, yields 



-,^ {Fy + 2TzX) +(r2^y2)-2 



(28) 



which possesses a local maximum at a; = FyTz/{2T^), 
y^ = 1 — a;^ (the sign of y being irrelevant). This is valid 
only for \FyTz\ < 2T^, otherwise x = sgn(f;f'^), y = 
has to be chosen. In the latter case, the value taken on 
is 



a 



K?-F4^ 



rp2 



(29) 



For the general case (occurring naturally in simu- 
lations), we resorted to locating the maximum of ((25)) 
numerically (by means of a Lagrangian multiplier, the 
Newton-Raphson method and being careful about poles) . 
Let us term this and the subsequent usage of the max- 
imized <P in (1^^ the application of the full (Tresca- 
based) criterion. In section 31 we will come back to it 
and compare the results to the ones obtained by em- 
ploying the following rather crude estimate. 

To obtain an approximate closed formula, we disre- 



gard the exact constraint x 



1 and replace it 



hy x^ — y^ = \. Then again, the lack of local maxima 
in enforces its maximum to be located at one of the 
four corner points {x, y) = (±1, ±1). This in turn, given 
the form of (pSj). shows that the signs of x and y have 
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to be those of Fy/T^ and F^T^ — 2Fj;Tz, respectively. 
Plugging that back into ([M)) yields 



F 



Trcsca 



^ F^ 



{2\F,T, 



f F2/4 + 4(T2 + 2T2)/a2 
2F,T,\+4\FyT,\)/a . 



(30) 



Comparing it to the special cases (051) . (071) . and (1^ . 
we recognize that the estimation can be improved by 
just dropping the factor 2 in front of T^, i.e. we have 
to check 



F 



Tic 



4 



.T2 + r2 



a 



(31) 



against the threshold F^. Let us call this the simplified 
criterion. 

For comparison purposes, we introduce the decou- 
pled criterion as well, where all loading modes are checked 
individually: 



■\^F^ + Fl |F,|/2, 2\T^\/a, 2\T^\/a} 



(32) 

Whatever criterion (i.e. definition of Fxresca) used, 
in the case of exceeding the threshold (i.e. -Frresca > 
-F*), the contact is irreversibly transformed into a usual, 
frictional linear spring-dash pot one with the same visco- 
elastic parameters. 



3 Simulation 

3.1 Simulation setup 

In this part, we describe in detail the procedure which 
was used to create standard yield tests on cylindrical 
shaped specimens using the LAMMPS [10 molecular 
dynamics simulation code, extended to treat caked con- 
tacts according to the model in section [2] 

In order to create a homogeneous, random close 
packing, first particles were poured into a half cylinder 
placed horizontally in a vertical gravity field. The cylin- 
der had the same diameter as the desired test specimen 
but was necessarily much longer. 

Particles had a narrow size distribution that is large 
enough to prohibit crystallization on the cylinder wall. 
The number of particles varied from 200 to 10000, the 
friction coefficient in the preparation as well as in the 
tests was set to 0. 

The radius of the cylinder was chosen to be 



R 



' N 



1/3 



(33) 






(a) 



(b) 



(c) 



Fig. 1 (a) A typical test setup witli aspect ratio s 3 and 
N = 200. Tlie smootli, blue particles are the driving particles, 
the striped brown ones are the bulk particles, (b) The contact 
network of the same system, (c) The broken contact network 
after strain e = 0.00017. Note, that particle displacement 
would not be visible. 



where d is the average particle diameter, N is the num- 
ber of particles and s is the desired height width sam- 
ple ratio: s=H/{2R), for which values in the range of 
0.5—10 were chosen. The value of / describes the aver- 
age volume taken up by a particle as 2d^ j f and thus, 
for a random close packing with a volume fraction of 
around 0.65, amounts to / ~ 2.5. 

After all particles were poured into the cylinder the 
sample was uni-axially compressed by two plates in ab- 
sence of gravity. The compression force employed was 
small enough to keep the average deformations of the 
friction-less, linear spring-dash pot contacts below 10~^ 
particle radii. For the subsequent failure tests, the stiff- 
ness was reduced by one order of magnitude and the 
pairwise distance was rescaled to avoid pre-stressing. 

The tests were done in absence of gravity and with- 
out the compressing plates. There are fewer contacts 
from the bulk to a plate than to a layer of particles, thus 
if plates were used for driving it would be expected that 
the contacts between a plate and the touching particles 
would fail first. So we chose the last layers of particles 
(the ones that are closer than 1.1 particle radius to the 
compressing plate) to move as a rigid body and drive 
the system. This allows the failure zone to be every- 
where in the sample. A typical setup and test result is 
shown on Fig. [1] 

Two different driving methods were used: (i) the 
driving layer had infinite mass, (ii) the driving layer 
had infinite mass only in the driving direction but fi- 
nite in the perpendicular directions. Since no significant 
difference was observed, the results reported in this pa- 
per are shown with driving particles of infinite mass. 
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The driving velocity was always kept low enough as to 
avoid kinetic effects. 

After the preparation the sample was tested. At the 
beginning of the tests all existing (i.e. force transmit- 
ting) contacts were "glued" , all subsequent contacts (in- 
cluding broken ones at later times) were handled with 
the standard Hooke contact law including static fric- 
tion. The simulation was run long enough such that 
the upper and lower part are completely separated as 
shown on Fig. [ijc). 

The different tests were realized by giving the driv- 
ing particles on both sides of the sample initially a ve- 
locity (or angular velocity) in the required direction but 
with opposing sign. 



3.2 Parameters 

In order to interpret the results of numerical simula- 
tions we have to take a closer look at the dimensions 
of the parameters. Due to the quasi static nature of 
the tests (up to the failure point) and the absence of 
gravity, quantities like time and mass are unimportant, 
instead force and distance units determine the scales of 
the numerical results. 

Two parameters of dimension distance enter: the av- 
erage particle diameter d and the contact glue diameter 
2a. Since in all simulations the ratio a=0.1d was used, 
d is taken as the unit of length. 

The stress or force scale is set by the stress thresh- 
old cr* (cf. (EH)) and i^, — na'^at,, respectively. Due to 
technical reasons, 17 times the latter was chosen as unit 
of force. 

Thus, from now on all quantities will be given in 
this natural units {d for length and 17i^* for force). 

The normal contact stiffness fc„ (the tangential one 
is set to kt=kn/2, corresponding to the limit = 0, 
cf . (IT4l) ) relates force and length and must fulfill fc„ ^ 1 
to assure the brittle limit. This can also be expressed 
as the failure strain scale 



krt.d 



(34) 



being very small. 

The preparation was done with the aim to achieve 
a sample as homogeneous and isotropic as possible. For 
such a random configuration, the macroscopic Young's 
modulus can be calculated analytically [5J[7] from the 
stiffness parameters of the contacts: 



E : 



Nj{2kn+3kt) 

W 4fc„ + kt 



(35) 



where Nc is the number of contacts and V is the volume 
of the sample. The latter can be calculated using (1551) : 

^ ^ z/(2fc„ + 3fc,) 

12d(4fc„ + kt)' ^ ' 

where z is the average coordination number which was 
measured to be z~5.06 in all cases. This value is less 
than the isostatic limit for homogeneous packing due 
to the boundary effects. Using the values /=2.5 and 
kt—kn/2, we get an estimate for the macroscopic Young's 
modulus: 



E = 0.81- 



(37) 



The macroscopic Poisson ratio can also be expressed as 
function of the contact stiffness [7J[5]: 

4fc„ + kt ^ ^ 

The actual simulation parameters (in natural units) 
chosen for the current study were taken to allow reli- 
able and fast simulation of the tests and yet fulfill the 
requirements above: 



d= 1 
fc„ = 1256 
a* = 1.88 



a = 0.1 



(39) 



From those we can calculate the macroscopic Young's 
modulus and Poisson ratio ([55]) as well as the limit 
strain (l34l): 



E = 1017 



I' ^ 0.11 



9.4 X 10"^ 



(40) 



4 Yield tests 

4.1 Poisson ratio 

The easiest quantity that we can measure in our sim- 
ulations is the Poisson ratio. In pull experiments we 
recorded the displacement of the driving particles in the 
z direction and the radius of the sample at z^H/2. The 
latter was done by averaging the outer radius over all 
surface particles that have at least two contacts (remove 
rattling, dangling particles) . Since the relative displace- 
ment is very small (of the order of 10^^) we can use the 
linearized expression for the Poisson ratio: 

..-If (.1) 

We averaged for different samples and system size to 
get the result: 



ly = 0.108 ± 0.01 



(42) 



This value is in perfect agreement with the prediction 
of (gOl). 
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4.2 Stress/strain behavior 

Three different tests were done: pull, shear and twist. In 
the tests the driving was realized by giving the driving 
particles a constant (angular) velocity. The net force 
(torque), the number of intact contacts, the deforma- 
tion of the sample was recorded as function of the dis- 
placement of the driving particles. 

The stress/strain evolution of the samples in the 
homogeneous and isotropic case are described by three 
parameters: Young's moduluf^ E, Poisson ratio v and 
maximal strain e* . So the evolution of all test should be 
described by one single curve. Thus the measured dis- 
placement and force (torque) data has to be converted 
to a stress/strain relation. 

It is more informative if this conversion is done with 
the dimensionless quantities of the aspect ratio s and 
the number of particles N since it allows immediately 
the scaling of different tests in a transparent manner. 

In general the stress and the strain are not homo- 
geneous, so we focus on their two important character- 
istics: the maximum and the average. The maximum 
stress governs the point where the first contact breaks, 
the average value indicates the importance of the in- 
homogeneities of the stress and we compare it to the 
maxima of the stress/strain curve. 

The pull deformation is homogeneous thus the max- 
imum is the same as the average: 



Z 



(/vr)V3z 



£p-ep,max- ^ - 2s2/3^rl/3j' 



(43) 



where Z is the vertical displacement of the driving par- 
ticles and / is the volume factor of ([55]) . The stress can 
be taken directly from section 12.11 



2/3 



El 



(44) 



The twist deformation has a radius dependence, but 
apart from this is very similar to the previous case. 

Ra a 
H Is 

et = 2et,,nax/3, (45) 

where a is the angle covered by the driving layer. The 
form of the stress again is identical to the microscopic 
one in 12.11 

_ 4r(l + z^) _ 4/(l + zy)s T 



o-t,, 



(ft = 2crt^niax/3 



N 



J3 



(46) 



^ Just like in the previous section, we use the symbols E and 
u for the macroscopic quantities instead of the ones at contact 
level as in section [2] The same apphes to the deformations 
X, Z, to forces/torque F^, F^, T and stresses. 



The shear deformation is a bit more involved than 
the previous ones, since its behavior is fundamentally 
different whether the aspect ratio is large or small com- 
pared to unity (with a crossover in between). The short 
cylinder limit was covered in Sec. 12.11 The long cylin- 
der limit is more complicated, but consists of a combi- 
nation of bending and shearing[T3] as well. We define 
here the strain as the most straightforward interpola- 
tion of the two cases that recovers the short cylinder 
result for s — ?> and the long cylinder form in the limit 
s — > oo: 



X 



4£t 



Stt 



(47) 



Applying Hooke's law to the whole strain field yields 
the interpolating stress 



2zy + 2 



3W 



7ri?2(l 



fs 



2/3 (2u + 2 + 



(1 -I- si A) d? 



4cr., 



Stt 



(48) 



Figured] (a) shows the average stress/strain curves 
for different experiments, system size, sample length. 
The measured Young's modulus varies around a mean 
value _E~880, which is 15% less than the theoretical 
value of ([37| . There is a slight scatter of the linear slope 
of different systems but in spite of the moderate system 
size still one can observe a convincing data collapse with 
different system size, aspect ratio and test type. 

In Fig.[3](b) the maxima of the average stress/strain 
curves are shown as function of the average strain. It 
is evident that the elastic behavior is very close to the 
Hooke's law and can be described by a single Young's 
modulus. The average macroscopic stress and strain are 
always smaller than the scale (l40l) , as must be expected 
due to spatial disorder. 

The maximal stress/strain values at the moment 
of the first broken contact are plotted in Fig. [3] (c). 
The values are scattered around e* of (|i(71). We note 
two important features: (i) points which are above 
are mainly for system size iV=200 and for shear and 
twist tests. This is the result of the inhomogeneities of 
our system. Due to the preparation procedure the vol- 
ume fraction at the lateral surface of the cylinder is 
larger than in the bulk. This was tested by repeating 
a few tests with a sample cut out from a bigger one, 
and this feature disappeared, (ii) Maximal stresses for 
long cylinders in shear tests are somewhat larger than 
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0.16 



Fig. 2 Color online, (a) Stress strain evolution (b) The max- 
imum of the average stress, (c) End of elasticity (first broken 
contact) in the stress/strain plain. Colors represent experi- 
ment type pull: blue, shear: green, twist: red. Symbol fill in- 
tensity indicates system size N: 200; half, 2000: empty, 10000: 
full. Symbol types show aspect ratio s: 1: o, 2: □, 3: Q), 5: A, 
7: V, 8: +, 10: x 



expected, these are the resuh of the same skin effect 
during the bending of the sample. 



4.3 Yield criteria 

In the simulations we used all three criterions intro- 
duced in Sec. 12.21 We performed test using exactly the 




Fig. 3 The stress and the strain points where the first con- 
tact was broken. Colors represent experiment type pull: blue, 
shear: green, twist: red. Symbols indicate the different break 
criterions: Full -|-, Simplified: x. Decoupled □. System size 
iV=200, 2000, 10000, are from bottom to top respectively 
with artificial shifting of 0.3. The aspect ratio was chosen to 
be s=2 



same initial configuration with the different break cri- 
terions. Fig. 13] shows the summary of the results. 

The results show hardly any difference between the 
two coupled criterions (simplified (I5T]) and full (I^SI) . 
however the decoupled criterion allows the system to 
evolve without failure beyond the point where the oth- 
ers already broke. 

It is remarkable that the simple coupling of the 
different deformation modes leads to indistinguishable 
failure points for large systems. For small systems the 
only difference is in the twist deformation. 

The elapsed time during simulations was also mea- 
sured. We observed no measurable difference in simu- 
lation duration of tests with decoupled and simplified 
criterions. Simulations with the full criterion on average 
require 55% more time independently of system size. 
This advocates the use of the simplified criterion if com- 
putational time is critical but this overhead is relatively 
small in spite of the numerical maximization of (j25p 
which is in general affordable if precision is required. 



4.4 Failure process 

In this section we study the process of the failure in the 
tests. Our simulations allow the study of each individual 
broken contact. The question we want to answer here is: 
Which is the deformation type that is the most relevant 
for a give test. 

We identify 6 different deformation types which are 
the core of the simplified break criterion (l5lT) : 



tensile 
shear 
tilt 



F2/4 
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Fig. 4 Color online. Average (sample and time) deformation forces in just broken contact bridges scaled by the maximal force 
versus Ni,. (a) Pull test with s=3 (b) Twist test with s=3 (c) Shear test with s=3, (d) Shear test with s=0.5. The symbols 
represent the following deformation component: tensile: +, shear: x, tilt: •, torsion: □, tilt+tensile: ■, torsion+shear: Q- On 
the right hand side of the plots a typical example of force network plot is shown for the given system after the test. 



torsion : 4T|/a^ 
tilt + tensile 2\F;,t^ - 2F^T,ya 
torsion + shear 4\FyT^\/a, (49) 

In order to study the relevance of the above defor- 
mation types, the following procedure was used: At the 
moment of a contact break, the above quantities are 
recorded. 

Since we have a random packing, the values show 
huge variations from contact to contact. An average 
is necessary to be able to observe the different trends. 
We averaged the recorded Umit forces for each succes- 
sive 5% of broken contacts. So if Nt,{t) is the number 
of broken contacts at time t then the x axes shows 
nb=Nb{t)/Nb{T), with T being the time at the end of 
the experiment. This representation allows us to aver- 
age also for different sample realizations (5-10) so that 
early and late stages of failure can be identified, as 
shown in Fig. 

The beginning of the pull test [Fig. H] (a)] is marked 
by contacts strained in the tensile direction. This dom- 
inance is clear for the first 25% broken contacts but re- 
mains the most important up to 60% broken contacts. 
The end stage (after 75% of the contacts were broken) , 
is surprisingly characterized by contacts breaking due 
to torsion. Our view is that this is due to breaking of 



small chains, namely a single particle having two con- 
tacts one to each side of the sample. In general, since all 
"strong" contacts are already broken, its two contacts 
are loaded mainly with a combination of torsion and 
tilting. A remarkable feature is that tilting gets stabi- 
lized by shear and tensile forces, while torsion is not. 
This can be seen by the fact that tilting-t-tensile defor- 
mation is almost the second most important loading at 
late stages. Since we have used a rather narrow contact 
diameter a=0.1d, contacts break easily under torsional 
load. 

Twist tests [Fig. H (b)] show a surprisingly similar 
picture to pull experiment, except that at early stages 
it is the shear deformation which dominates. The late 
stage is also dominated by torsional deformation. 

The behavior of shear tests are fundamentally dif- 
ferent depending on the aspect ratio. If the aspect ratio 
is small [s < 1, Fig. 2] (d)] the overall behavior is sim- 
ilar to the twist deformation: shear dominates at early 
stages and torsion at late ones. 

It is common to all of the above cases that the 
stress/strain behavior of the sample is determined by 
the contacts that broke in the early stages, moreover 
the whole fracture is realized very fast in a cascade like 
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Fig. 5 (Color online) The number of contacts as function 
of the average strain for s=3 and A''=2000, for the tests: pull: 
blue, shear: green, twist: red 



manner as shown on Fig. [5j as to be expected in the 
brittle Hmit. 

The only test where a different scenario can be ob- 
served is the shear test with long cylinders [Fig. |4] (c)]. 
At earfy stages the tensile deformation dominates. This 
is due to the bending of the sample. The force network 
shows that the contacts start to break at two points si- 
multaneously. These points are the ones where the ma- 
terial is stretched the most (right at the moving "front" 
of the driving layer) . The stress is very localized so de- 
spite the damage the sample has still a significant re- 
sistance which is manifested by the long oscillating tail 
of the stress/strain curve [Fig. [2] (a)] and the slow step- 
wise decreasing of the number of contacts. 



decoupling the loads leads to a rather large error of 
25%, whereas the simplified criterion agreed with the 
full one within statistical errors. This applies to our 
tests, but in the case of the full criterion turning out to 
be necessary in other situations, the found overhead of 
55% for the numerical maximization is affordable. 

As a benchmark, we have done pull, shear and tor- 
sion tests with cylindrical specimens created from glued 
particles. These exhibited, up to the failure point, an 
elastic behavior with the same Young's modulus for all 
experiments and test parameters. This enabled us to 
collapse all different test results onto one single curve. 
The deviations of the theoretical prediction for the spec- 
imens' elastic properties were below 15%. 

During the tests we payed special attention to the 
dominant loading of the contacts that are about to 
break. We found that the first half of broken contacts 
in pull tests mostly break due to tensile stress, in twist 
tests due to shear stress, and in shear tests of very short 
samples due to shear stress as well. Long samples be- 
have differently in shear tests, though, due to bending 
the dominating load gets tensile instead of shear. Fi- 
nally, the late stages of the failure is mainly character- 
ized by contacts breaking under torsional load due to 
break up of twisted weak chains connecting the almost 
separated parts. 

For the sake of a clearer study of this relative im- 
portance of the loading types and the different failure 
criteria, we used a constant diameter of the contact 
cylinders. In future works, the stochastic nature of the 
geometry and strength of the contact, possibly obeying 
a wide distributioujGj, will be taken into account. 



5 Conclusion 



In this paper we introduced a new model for DEM sim- 
ulations with breakable cemented contacts. Contacts 
were considered as small, flat, elastic cylinders which 
keep the particles in contact until the shear stress in 
the cementing cylinder reaches a maximal value con- 
forming to the Tresca criterion of failure. This allowed 
for deriving the elastic and the failure behavior of the 
contact on a common basis. In the limit of very flat 
cylinders, the stiffness matrix gets decoupled, jus- 
tifying this common a priori modelling Ansatz. 

The problem of flnding the maximal shear stress 
in an arbitrarily loaded cylinder turned out to possess 
analytic solutions too complicated for practical pur- 
poses (but is numerically well accessible). Therefore 
we considered two approximations: The decoupled crite- 
rion coincides with the exact maximization only for sin- 
gle loading modes, while the simplified criterion agrees 
with the latter in more general situations, due to tak- 
ing into account load combinations. We showed that 
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